AD-A229  671 


one  FILE  COHY  ^ 

NASA  Contractor  Report  187460 
ICASE  Report  90-76 


ICASE 


TOWARD  THE  LARGE-EDDY  SIMULATION  OF 
COMPRESSIBLE  TURBULENT  FLOWS 


G.  Erlebacher 
M.  Y.  Hussaini 
C.  G.  Speziale 
T.  A.  Zang 


Contract  No.  NAS  1-18605 
October  1990 


Institute  for  Computer  Applications  in  Science  and  Engineering 
NASA  Langley  Research  Center 
Hampton,  Virginia  23665-5225 

Operated  by  the  Universities  Space  Research  Association 


fW\SA 

National  Aeronautics  and 
Space  Administration 

Langley  Research  Center 

Hampton,  Virginia  23665-5225 


D  fSTHTBUnON  STATEIlCa<T 

A^i'ntoved  for  pobllo  ^ 

]>L»V>lbutioD  UaUflalMj  ,  ' 


A 


V) 

i- 


TOWARD  THE  LARGE-EDDY  SIMULATION  OF 
COMPRESSIBLE  TURBULENT  FLOWS* 


G.  Erlebacher^ 

ICASE,  NASA  Langley  Research  Center 

M.  Y.  Hussaini^ 

ICASE,  NASA  Langley  Research  Center 

C.  G,  Speziale^ 

ICASE,  NASA  Larigley  Research  Center 


T.  A.  Zang 

NASA  Langley  Research  Center 
Hampton,  VA  23665 


ABSTRACT 


'  New  subgrid-scale  models  for  the  large-eddy  simulation  of  compressible  turbulent  flows  are 
developed  and  tested  based  on  the  Favrc-filtered  equations  of  motion  for  an  ideal  gas.  A 
compressible  generalization  of  the  linear  combination  of  the  Smagorinsky  model  and  scale- 
similarity  model,  in  terms  of  Favre-filtered  fields,  is  obtained  for  the  subgrid-scale  stress  ten¬ 
sor.  An  analogous  thermal  linear  combination  model  is  also  developed  for  the  subgrid-scale 
heat  flux  vector.  The  two  dimensionless  constants  associated  with  these  subgrid-scale  models 
are  obtained  by  correlating  with  the  results  of  direct  numerical  simulations  of  compressible 
isotropic  turbulence  performed  on  a  96^ .'grid  using  Fourier  collocation  methods.  Extensive 
comparisons  between  the  direct  and  modeled  subgrid-scalc  fields  are  provided  in  order  to  val¬ 
idate  the  models.  A  large-eddy  sy^lation  of  the  decay  of  compressible  isotropic  turbulence 
-  conducted  on  a  coarse  32^^  grid  •*  is  shown  to  yield  results  that  are  in  excellent  agreement 
with  the  fine  grid  direct  simulation.  Future  applications  of  these  compressible  subgrid-scale 
models  to  the  large-eddy  simulation  of  more  complex  supersonic  flows  are  discussed  briefly. 
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1.  Introduction 


The  direct  numerical  simulation  of  turbulent  flows  at  the  high  Reynolds  numbers  encountered 
in  problems  of  technological  importance  is  all  but  impossible  as  a  result  of  the  wide  range 
of  scales  that  are  present.  Consequently,  the  solutions  to  such  problems  must  invariably 
be  based  on  some  form  of  turbulence  modeling.  Traditional  turbulence  models  based  on 
Reynolds  averages  have  had  only  limited  success  since  the  large  scales  of  the  turbulence 
-  which  contain  most  of  the  energy  -  are  highly  dependent  on  the  geometry  of  the  flow 
being  considered.  Experience  has  indicated  that  such  models  usually  break  down  when  a 
variety  of  turbulent  flows  are  considered  (Lumley  1983).  The  small  scales  are  more  universal 
in  character,  and  serve  mainly  as  a  source  for  dissipation.  Hence,  it  can  be  argued  that  a 
better  understanding  of  turbulent  flows  could  be  achieved  if  just  the  small  scales  are  modeled 
while  the  large  scales  are  calculated  (Deardorff  1970).  This  is  the  fundamental  idea  behind 
large-eddy  simulations. 


During  the  past  decade,  considerable  progress  has  been  made  in  the  large-eddy  simulation 
of  incompressible  turbulent  flows.  This  effort  has  shed  new  light  on  the  physics  of  turbulence. 
The  earliest  work  relied  heavily  on  the  use  of  the  Reynolds  averaging  assumption  to  elim¬ 
inate  the  Leonard  and  cross  stresses  while  the  Reynolds  stresses  were  computed  using  the 
Smagorinsky  model  (Deardorff  1970,  Leonard  1974,  Reynolds  1976).  More  recent  large-eddy 
simulations  have  been  based  on  the  direct  calculation  of  the  Leonard  stresses  with  models 
provided  for  the  cross  and  Reynolds  subgrid-scale  stresses  in  order  to  enhance  the  numerical 
accuracy  (see  Biringen  and  Reynolds  1981,  Bardina  Ferziger  and  Reynolds  1983).  However, 
among  these  newer  models,  only  the  Bardina,  Ferziger  and  Reynolds  (1983)  model,  with  a 
Bardina  constant  of  1.0,  satisfies  the  important  physical  constraint  of  Galilean  invar’ance 
(Speziale  1985).  The  underlying  physical  concepts,  fundamental  numerical  algorithms,  and 
comprehensive  historical  data  behind  the  recent  field  of  large-eddy  simulation  have  been 
presented  in  articles  by  Schumann  (1975),  Yoke  and  Collins  (1983)  and  Rogalio  and  Moin 
(1984).  More  recently,  work  on  the  subgrid-scale  modeling  of  transition  to  turbulence  of  ini¬ 
tially  laminar  incompressible  flows  has  begun  (Piomelli,  Zang,  Speziale  and  Hussaini  1990). 
Several  large-eddy  simulations  have  been  performed  and  initial  results  are  promising. 


Despite  the  intensive  research  effort  that  has  been  devoted  to  the  large-eddy  simulation 
of  incompressible  flows  as  outlined  above,  it  appears  that  no  large-eddy  simulation  of  a 
compressible  turbulent  flow  has  yet  been  attempted.  Of  course,  such  work  could  have  im¬ 
portant  technological  applications  in  the  analysis  of  turbulent  supersonic  flows,  where  shock' 
waves  are  generated,  and  in  turbulent  flows  within  combustion  chambers.  The  prerequisite 
for  carrying  out  such  computations  is  the  development  of  suitable  subgrid-scale  models  for 
compressible  turbulent  flows.  With  the  exception  of  the  recent  work  of  Yoshizawa  (1986) 
and  Speziale  et  al.  (1988),  few,  if  any,  studies  along  these  lines  appear  to  have  been  pub¬ 
lished.  The  subgrid-scale  models  of  Yoshizawa  are  only  suitable  for  slightly  compressible 
turbulent  flows  since  they  made  use  of  an  asymptotic  expansion  about  an  incompressible 
state.  Recently  however,  Dahlburg,  Zang  and  Dahlburg  (1990)  have  performed  an  extensive 
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parameter  study  using  the  model  developed  by  Speziale  et  al.  (1988).  It  was  shown  by 
Speziale  et  al.  (1988)  that,  for  the  purposes  of  accuracy,  the  Leonard  and  cross  stresses 
must  be  accounted  for.  Furthermore,  the  modeling  of  the  isotropic  part  of  the  Reynolds 
subgrid-scale  stress  tensor  was  shown  to  be  questionable  -  an  issue  that  was  left  for  future 
research. 

In  this  paper,  complete  subgrid-scale  models  are  developed  for  the  closure  of  the  Favre- 
filtered  Navier-Stokes  and  energy  equations.  The  compressible  subgrid-scale  stress  model 
that  is  obtained  in  Section  2  reduces  to  the  linear  combination  model  of  Bardina,  Ferziger 
and  Reynolds  (1983)  in  the  incompressible  limit.  Likewise,  the  subgrid-scale  heat  flux  model 
that  is  obtained  herein  consists  of  an  analogous  linear  combination  of  scale  similarity  and 
gradient  transport  terms.  The  dimensionless  constant  which  appears  in  the  subgrid-scale 
stress  model  is  arrived  at  through  correlation  analysis  of  data  generated  from  direct  numerical 
simulations  of  compressible  isotropic  turbulence.  A  more  detailed  comparison  of  computed 
and  modeled  subgrid-scale  fields  is  presented  along  with  the  results  of  a  large-eddy  simulation 
of  compressible  isotropic  turbulence. 

2.  Subgrid-Scale  Models  for  Compressible  Turbulence 


The  compressible  turbulent  flow  of  an  ideal  gas  is  considered.  Such  flows  are  governed  by 
the  continuity,  momentum  and  energy  equations  which  -  neglecting  body  forces  -  are  give.i 
by  (cf.  Batchelor  1967) 


dp  d{pvk) 
dt  dxk 

(i) 
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respectively,  where  p  is  the  mass  density,  v  is  the  velocity  vector,  p  is  the  thermodynamic 
pressure,  pi  is  the  dynamic  viscosity,  h  is  the  enthalpy,  T  is  the  absolute  temperature,  and  k 
is  the  thermal  conductivity.  The  viscous  stress  cr*(  and  the  viscous  dissipation  are  defined 
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respectively.  Herein,  the  Einstein  summation  convention  applies  to  repeated  indices.  Equa¬ 
tions  (l)-(3)  must  be  supplemented  with  the  equations  of  state 


P  -  pRT,  h  -  CpT 


(6) 
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for  an  ideal  gas  where  R  is  the  ideal  gas  constant  and  Cp  is  the  specific  heat  at  constant 
pressure.  Likewise,  the  dependence  of  the  viscosity  and  thermal  conductivity  on  the  temper¬ 
ature  must  be  provided  (i.e.,  relationships  of  the  form  fi  =  /r(T)  and  k,  =  x^iT)  are  needed 
and  these  depend  on  the  gas  under  consideration). 

Any  flow  variable  :F  can  be  filtered  in  the  following  manner: 

J(x)=  /  G(x-z,A)J^Wd^z  (7) 

where  G  is  a  filter  function,  A  is  the  computational  mesh  size,  and  D  is  the  domain  of  the 
fluid.  The  filter  function  G  is  typically  taken  to  be  an  infinitely  differentiable  function  of 
bounded  support  in  a  bounded  domain,  or  a  Gaussian  distribution  in  a  periodic  domain^  It 
is  normalized  by  requiring  that 

/  G(x-2,A)d^2  -  1.  (8) 

JD 

It  follows  that  in  the  limit  as  the  computational  mesh  size  goes  to  zero,  (7)  becomes  a  Dirac 
delta  sequence,  i.e. 

hm^y  G(x  -  z,  A)J^(z)(i^z  =  J  6{-k  -  z)f{z)d^z  —  (9) 

where  ^(x-z)  is  the  Dirac  delta  function  (Arfken  1970).  The  filter  function  has  the  property 
that  the  amplitude  of  the  high-frequency  spatial  Fourier  components  of  any  flow  variable  ^ 
are  substantially  reduced.  Consequently,  represents  the  large-scale  part  of  !F.  At  this 
point,  it  should  be  mentioned  that  as  a  result  of  the  defining  properties  of  G,  it  follows  that 

dt  dt  ’  dxk  dxk 

Piomelli,  Ferziger  and  Moin  (1987)  discuss  the  relationship  between  the  form  of  the  filter 
function  and  that  of  the  subgrid-scale  turbulence  model. 


The  turbulent  fields  are  decomposed  as  follows,  based  on  Favre  filtering: 

=  (11) 
where  the  Favre  filter  _ 

(12) 

P 

is  defined  in  an  analogous  manner  to  the  Favre  time  average  which  has  been  of  use  in  the 
more  traditional  studies  of  compressible  turbulent  flows  (Hinze  1975).  However,  contrary  to 
the  more  traditional  Favre  time  averaging. 


T 

*A  Gaussian  filter  is  adopted  for  the  calculations  in  this  study 
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in  general,  and  hence 


(14) 


0- 

The  direct  filtering  of  the  continuity  equation  (1)  yields 

dp  ,  d[pvk) 


dt  dxk 


=  0 
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where  we  have  used  (10)  and  (12).  Likewise,  a  direct  filtering  of  the  momentum  equation 
yields 

d{pvk)  d{pVkVi)  dp  ,  daki  ,  dr^ 
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where 

and 
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Tkl  =  -p{VkVl  -  VkVl  +  WfcV,  +  v[Vk  +  v'^v\)  (18) 

is  the  subgrid-scale  stress  tensor.  The  subgrid-scale  stress  tensor  can  be  decomposed  as 
fellows. 


where 


r  =  L-i-C  +  R 

(19) 

Lki 

=  -p{VkVt  -  V*V/) 

(20) 

Cki 

=  -p{v'^vi  -t-  v\vk) 

(21) 

Rki 

=  -pv'kV'l 

(22) 

are  respectively,  the  subgrid-scale  Leonard,  cross,  and  Reynolds  stresses  based  on  P'avre 
filtering.  From  (20),  it  is  clear  that  the  Leonard  stress  can  be  calculated  directly  and  does 
not  need  to  be  modeled.  The  cross  stress  is  modeled  with  the  scale  similarity  model 


Ckt  =  -fivkVt  -  V^V[) 


(23) 


(with  a  cucfficient  of  unity  to  ensure  Galilean  invariance  of  the  overall  model).  This  model 
is  analogous  to  its  incompressible  counterpart,  which  has  been  reasonably  successful  in  the 
large-eddy  simulation  of  incompressible  turbulent  flows  (Bardina,  Ferziger  and  Reynolds 
1983,  Speziale  1985).  The  subgrid-scale  Reynolds  stress  tensor  is  separated  into  deviatoric 
and  isotropic  parts,  respectively,  as  follows: 


where 


and 


R  =  dR  +  /R 


D^ki  =  -pWkv'i'-  -^vy^Ski), 


iRkt  =  -■^p'^yxhi- 
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Here,  the  deviatoric  part  of  the  subgrid-scale  Reynolds  stress  tensor,  oR-i  is  modeled  using 
the  compressible  generalization  of  the  Smagorinsky  model  that  is  given  by 


nRki  =  2CRpA^Iiy^{Skt  —  --SmmSki) 


(27) 


where 


a  .  dvi^ 

“  2^dxi  dxj 


lie  =  SmnS, 


(28) 

(29) 


(i.e.,  S  is  the  Favre  filtered  rate  of  strain  tensor  while  II§  is  its  second  invariant)  and  Cr  is 
the  compressible  Smagorinsky  constant.  Yoshizawa  -  by  means  of  a  two-scale  DIA  method  - 
derived  a  model  for  the  isotropic  part  of  the  subgrid-scale  Reynolds  stress  tensor,  /R,  given 
by 

iRkt  =  -lCipA^Il5S„  (30) 

where  C/  is  a  dimensionless  constant.  Equation  (30)  can,  for  the  most  part,  be  obtained 
by  making  a  turbulence  production  equals  dissipation  equilibrium  hypothesis  (Yoshizawa 
1986).  However,  this  model  was  shown  by  Speziale  et  al.  (1988)  to  correlate  very  poorly 
with  the  results  of  direct  numerical  simulations  of  compressible  isotropic  turbulence.  Since 
/R  is  extremely  small  compared  to  the  thermodynamic  pressure,  we  propose  to  neglect  it  - 
an  assumption  that  will  be  juslifed  later.  Hence,  the  overall  subgrid-scale  stress  model  we 
propose  takes  the  form 


Tki  =  -p{vk^i  -  WfcV/)  +  2CRpA'^Iiy^{Ski  - 


l-r, 


M- 


(31) 


In  the  incompressible  limit,  Eq.  (31)  reduces  to  the  linear  combination  model 


-  rulp  =  UfcUi  -  VkVi  -  2CflA*//y*5fci 


(32) 


of  Bardina,  Ferziger,  and  Reynolds  (1983)  v/here  the  Bardina  constant  is  one  in  order  to  sat¬ 
isfy  Galilean  invariance  (Speziale  1985).  This  reduction  process  is  a  consequence  of  Eq.  (31) 
and 

V  =  V,  Smm  =  =  0  (33) 

when  p  becomes  constant. 


A  direct  filtering  of  the  energy  equation  yields  the  filtered  form 


a(c,pr)  ,  d(c,pi,f)  dp  ,  Sf  ,  ^ 

d  dQ, 

dxk  dxfc  dxk 

where  __  _ _ 

Qk  =  Cj,p  {vkT  -  Vkf  +  v'lj'  +  VkT'  +  v'^T') 
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is  the  subgrid-scale  heat  flux.  The  subgrid-scale  heat  flux  can  be  decomposed  in  the  same 
fashion  as  the  subgrid-scale  stresses.  This  leads  to 


where 


Q  =  -t- 

(36) 

=  C^p{ii^-vf) 

(37) 

=  C^p(<?’ +  uTT') 

(38) 

=  CpPv'T' 

(39) 

are  the  Leonard,  cross,  and  Reynolds  heat  fluxes.  Analogous  to  the  modeling  of  the  cross 
stress,  the  cross  heat  flux  is  modeled  using  the  scale  similarity  format 

f).  (40) 

The  Reynolds  heat  flux  is  modeled  with  the  usual  gradient  transport  format  as  follows  (cf. 
Eidson  1985): 

(41) 

where  Prj  is  the  turbulent  Prandtl  number.  Of  course,  the  Leonard  heat  flux  can  be 
calculated  directly.  Hence,  the  overall  model  for  the  subgrid-scale  heat  flux  we  propose  is  as 
follows: 

=  c,p  i(vTf  -  i.f)  -  ^ .  («) 

and  is  obtained  by  combining  equations  (37),  (40),  and  (41). 


At  this  point,  some  comments  need  to  be  made  concerning  the  viscous  terms  on  the 
right-hand  side  of  (16)  and  the  pressure  gradient-velocity  and  viscous  dissipation  terms 
which  appear  on  the  right-hand  side  of  (34).  The  pressure  gradient-velocity  correlation  can 
be  written  in  the  alternative  form 

dp  _  d(pv]^)  ^  dvk 

^dxk  dxk  ^  dxk 
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where  only  the  pressure  dilatation  term 
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dxk 

model  and  not  much  success  has  been  achieved  in  dealing  with  it  in  the  context  of  Reynolds 
stress  models.  However,  within  the  framework  of  subgrid-scale  modeling,  this  term  and  its 
corresponding  cross  correlation  have  physical  interpretations.  They  represent  the  contribu¬ 
tion  of  the  dilatation  of  the  small  scales  to  the  internal  energy  variation  of  the  fluid  -  an 
effect  which  is  expected  to  be  small.  Hence,  for  this  initial  study,  we  neglect  such  terms. 
Furthermore,  since  the  mean  temperature  is  constant  and  the  temperature  fluctuations  are 
small  ( <  10%),  the  viscosity  and  thermal  conductivity  are  held  constant.  For  similar  reasons, 
we  also  neglect  the  small  scale  component  of  the  viscous  dissipation. 

The  turbulence  model  proposed  herein  is  thus  complete  once  values  for  the  constants 
Cr  and  Prx  are  obtained.  This  will  be  accomplished  using  the  results  of  direct  numerical 
simulations  of  compressible  isotropic  turbulence. 


)  is  extremely  difficult  to 


is  not  yet  closed.  The  temperature  dilatation  correlation  (T' 


3.  Numerical  Method 


Our  direct  simulations  of  compressible  turbulence  are  based  on  a  non-dimensional  form  of 
Eqs.  (l)-(3),  with  the  time  derivative  in  the  energy  equation  written  solely  in  terms  of  the 
pressure.  In  order  to  alleviate  the  severe  stability  limit  imposed  at  very  low  Mach  numbers 
by  the  acoustic  waves,  a  splitting  method  is  adopted.  The  first  step  integrates  the  equations 


(45) 


djpvk)  dipvkvi)  _  doki 

dt  dxi  dxi  ’ 
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—  +  7Pa 
dt  dxk  dxk 
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+  (7  -  1)^, 
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while  the  second  step  integrates 

^  d{pvk)  ^ 
dt  dxk 
djpvk)  .  dp_  _ 
dt  dxk 


(48) 

(49) 


^  ,  2^(P^fc)  _ 
dt  dxk 

The  constants  Cq  and  M^o  are  the  current  root  mean  square  (rms)  value  of  the  sound  speed  (c) 
and  the  reference  Mach  number,  while  7  =  Cp/Cv,  where  is  the  specific  heat  at  constant 
volume.  These  equations  are  non-dimensionalized  in  terms  of  a  length  scale  {Lo),  a  velocity 
scale  {Uq),  a  pressure  scale  (Pq),  a  reference  viscosity  p  and  a  reference  thermal  conductivity 
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K.  The  Reynolds  number  is  given  by  Re  —  pqUoLqIh,  the  Prandtl  number  by  Pr  —  Cppln, 
and  the  reference  Mach  number  by  Moo  =  Uoli^RToY^^ .  For  all  calculations  presented  in 
this  study,  7  =  1.4,  and  Pr  =  0.7.  Initially,  the  density  po  is  uniform  and  equal  to  one.  The 
computational  domain  is  a  cube,  normalized  to  (0,27r]^.  Periodic  boundary  conditions  are 
imposed  in  all  three  directions. 


The  spatial  derivatives  in  these  equations  are  approximated  by  a  Fourier  collocation 
method  (see,  for  example,  Hussaini  and  Zang  1987).  In  each  coordinate  direction,  N  grid 
points  are  used:  Xk  =  2Trj/N,  for  j  ~  0, 1, N  —  1.  The  derivative  of  a  function  P'{x)  with 
respect  to  Xk  is  approximated  by  the  analytic  derivative  of  the  trigonometric  interpolant  of 
.^(x)  in  the  direction  i*.  Most  simulations  of  incompressible,  homogeneous  turbulence  have 
used  a  Fourier  Galerkin  method.  The  compressible  equations,  however,  contain  cubic  rather 
than  quadratic  nonlinearities  and  true  Galerkin  methods  are  more  expensive  (compared  with 
collocation  methods)  than  they  are  for  incompressible  flow.  The  essential  difference  between 
collocation  and  Galerkin  methods  is  that  the  former  are  subject  to  both  truncation  and 
aliasing  errors,  whereas  the  latter  have  only  ti  n  errors.  As  discussed  extensively  by 

Canuto,  et  al  (1987),  the  aliasing  terms  are  n  "it  for  a  well-resolved  flow.  However, 

care  is  needed  to  pose  a  collocation  method  iri  ?  which  ensures  numerical  stability.  For 

this  reason,  the  second  term  in  (46)  is  actually  co  ;  the  equivalent  form 


d{pVkVi)  ,  dvk  d{pvi) 

— T +  pvi- h  V*-— 

UXi  (JXi  OXi 


(51) 


As  noted  by  Feiereisen,  Reynolds  and  Ferziger  (1981),  when  this  form  is  employed  together 
with  a  symmetric  differencing  method  in  space  (for  example  Fourier  collocation),  then  in 
addition  to  mass,  and  morrientum,  energy  is  conserved  for  the  ideal  compressible  equations 
(zero  viscosity  and  thermal  conductivity)  in  the  absence  of  time  differencing  (and  splitting) 
errors. 


The  second  fractional  step  of  the  splitting,  given  by  (48)-(50),  contains  most  of  the  effects 
of  the  acoustic  waves.  This  splitting  is  employed  at  each  stage  of  a  third-order  Runge- 
Kutta  method.  In  the  simulations  reported  here,  the  second  fractional  step  is  integrated 
analytically.  In  Fourier  space,  (48)-(50)  become 


^  -I-  ikirhi  =  0 


(52) 


dmi 

dt 


+  ikip  =  0 


(53) 


^  -f  iclkimi  =  0  (54) 

where  mi  =  pvi  and  Fourier  transformed  quantities  (which  depend  upon  the  wavenumber  k) 
are  denoted  by  a  circumflex.  The  exact  solution  of  these  equations  is 


—  \A  cos(coA:At,)  +  D  siii(coA:AtJ  -  A] 
^0 


(55) 


8 


(56) 


m 


W  - 


iki 


m|  ' - -{A  sin(coA:Afj)  —  B  cos{cokAt,)  +  B] 

Cok 

—  A  cos{cokAt,)  +  B  sin(coAt,) 


(57) 


where  k  -  |k|  is  the  magnitude  of  the  Fourier  wavenumber,  A  =  B  —  i^ki7h[^\  The 
superscript  1  denotes  the  result  of  the  first  fractional  step  of  the  splitting  and  the  superscript 
^  ♦he  results  of  the  second  fractional  step.  The  effective  time-ctep  of  the  Runge-Kutta  stage 
is  denoted  by  At,.  The  advantage  of  this  splitting  is  that  the  principal  terms  responsible  for 
the  acoustic  wa'^es  have  been  isolated.  Since  they  are  treated  semi-implicitly,  one  expects 
the  time-step  limitation  to  depend  upon  u  +  |c  —  co|  rather  than  v  +  c.  (Although  there  is 
also  a  viscous  stability  limit  for  the  first  fractional  step,  it  is  well  below  the  advection  limit 
in  the  cases  of  interest.)  This  is  clearly  a  substantial  advantage  at  low  Mach  numbers.  Since 
the  second  fractional  step  is  integrated  analytically,  it  does  not  contribute  to  any  time  step 
limitations.  If  one  is  truly  interested  in  all  the  details  arising  from  the  sound  waves,  or  if 
there  is  a  substantial  coupling  between  the  sound  waves  and  the  rest  of  the  flow,  then  the 
time-step  must  be  small  enough  to  resolve  the  temporal  evolution  of  these  waves.  But,  if 
only  the  larger-scale  sound  waves  are  of  interest,  then  this  splitting  method  is  useful. 


During  the  acoustic  fractional  step,  an  isotropic  truncation  is  performed:  for  each  variable 
{p,  p\  and  p),  all  Fourier  coefficients  for  which 

kiki  >  {NI2f  (58) 

are  set  to  zero.  This  reduces  the  numerical  anisotropy  produced  by  a  cubic  truncation. 
Moreover,  it  reduces  the  aliasing  interactions  in  the  collocation  method  (Canuto  et  al  1987, 
Chapters  3  and  7). 


The  compressible  code  can  also  be  executed  in  a  purely  explicit  mode.  In  this  case  no 
splitting  is  performed;  Eqs.  (45)-(50)  are  simply  combined  in  the  appropriate  manner  and 
integrated  directly. 


The  expected  stability  limit  of  this  three-dimensional  Fourier  collocation  method  for  the 
compressible  Navier-Stokes  equations  has  the  form 

■»  1  -J 


At  <  a 

where  for  the  semi-implicit  version, 


..j 


\N  J 


(59) 


U,  =  |u,|  d  |c  -  Col  (60) 

while 

-  |t'.|  +  |c|  (61) 

when  the  time  advancement  is  fully  explicit.  For  the  third  order  Runge-Kutta  method 
employed  here,  we  use  a  =  0.5. 


A  number  of  simulations  have  also  been  conducted  of  strictly  incompressible  flow.  These 
were  performed  v/ith  a  separate  code  which  also  used  a  Fourier  collocation  method,  but  for 
the  simpler,  incompressible  Navier-Stokes  equations. 
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4.  Comparison  with  Incompressible  Results 


The  initial  conditions  for  the  numerical  simulations  were  designed  to  reproduce  the  experi¬ 
mental  data  of  Comte- Bellot  and  Corrsin  (1971)  on  isotropic  turbulence,  hereafter  referred  to 
as  CBC.  These  experiments  were  also  the  basis  of  direct  simulations  used  by  Clark,  Ferziger 
and  Reynolds  (1979),  Bardina,  Ferziger  and  Reynolds  (1983),  and  McMillan  and  Ferziger 
(1979)  in  their  analyses  of  incompressible  LES  models.  Initial  conditions  are  chosen  to  match 
CBC  measurements  at  a  non-dimensional  time  of  240  (cf.  table  4  in  Comte-Bellot  and  Corrsin 
1971).  The  computational  domain  is  a  cube  of  side  20/27r  cm.  The  CBC  parameters  are 
associated  with  measurements  taken  behind  a  grid  with  a  mesh  spacing  of  one  inch,  and  a 
mean  fluid  velocity  of  393.7  in/sec.  The  initial  time  in  the  direct  simulation  corresponds  to 
t  =  0.00254  sec  in  the  CBC  experiment.  The  reference  length  {Lq),  velocity  (Uo)  and  pres¬ 
sure  (Po)  are  respectively  20/27r  cm,  1  cm! sec  and  1  gr/cmsec^.  After  generating  a  random, 
divergence-free  velocity  profile,  the  kinetic  energy  (in  Fourier  space)  is  scaled  to  match  the 
measured  CBC  energy  spectrum  (see  Appendix).  Finally,  the  velocities  are  scaled  so  that 
the  initial  rms  velocity  agrees  with  the  measured  values.  This  adjustment  is  typically  less 
than  1%,  which  provides  one  measure  of  the  uncertainty  in  the  fit  to  the  experimental  data. 
With  the  chosen  non-dimensionalization,  the  Reynolds  number  Re  =  UqLoIi^  is  22.74  based 
on  a  kinematic  viscosity  i/  =  0.14  crm? j sec.  Table  1  summarizes  the  parameters  measured 


CBC 

64^ 

96^ 

128^ 

With 

6.75 

6.75 

6.75 

6.75 

itr(Sk) 

- 

0.0 

0.0 

0.0 

E 

- 

68.3 

68.3 

68.3 

c 

462 

375 

432 

447 

■^11 

0.26 

0.28 

0.27 

0.25 

Ai2 

- 

0.20 

0.19 

0.27 

-^13 

- 

0.20 

0.19 

0.26 

Rx 

38.1 

43.6 

41.3 

37.8 

Table  1;  Initial  conditions  based  on  CBC  experiment  and  Clark  et  al. 
(1979)  calculation.  Mach  number  is  zero. 


by  CBC  at  t=240.  The  Taylor  microscale  length  is  defined  by 


^ki  - 


(62) 


!0 


and  the  dissipation  e  by 

c  =  2/i  J  S,jS,jd^x, 

where  5,j  is  the  rate  of  strain  tensor 


_l  (dv.  dvj 

2  [dx,  dx. 


(63) 


(64) 


In  Eq.  (62),  (•)  denotes  a  spatial  average.  Its  exact  definition  is  given  in  the  Appendix.  The 
Taylor  microscale  Reynolds  number  is 


Jix  = 


(65) 


The  velocity  derivative  skewness  and  flatness  tensors  Sk  and  FI  are  the  third  and  fourth 
moments  of  the  velocity  gradient  and  are  defined  by 


In  tables  1  and  2,  only  the  trace  of  the  sKewncss  and  flatness  tensors  are  shown.  The 
remaining  columns  list  the  parameters  obtained  from  the  initial  conditions  of  the  numerical 
simulations  on  64^,  96^  and  128^  grids.  There  is  a  20%  discrepancy  between  the  dissipation 
obtained  by  CBC  and  the  dissipation  computed  on  the  coarsest  grid  which  suggests  that  a  64^ 
grid  has  marginal  resolution,  at  best.  A  12%  difference  between  the  value  of  Rx  calculated 
on  the  64^  grid  and  that  obtained  by  CBC  confirms  the  need  for  grids  finer  than  64^.  On  a 
96^  grid,  both  the  dissipation  and  Rx  are  in  much  closer  agreement  with  CBC.  Discrepancies 
between  our  results  and  CBC  for  c  and  Rx  are  respectively  6.5%  and  7.5%  on  a  96^  grid.  On 
the  finest  grids  on  which  the  direct  simulations  were  performed,  the  computed  values  of  e 
and  Rx,  respectively,  have  relative  errors  of  3.5%  and  less  than  1%  when  compared  to  CBC. 


The  numerical  simulations  were  run  from  t  =  240  until  t  =  375  (in  CBC  units),  which 
corresponds  to  a  non-dimensional  tin"  interval  of  0.1145  (in  our  units).  Table  2  furnishes  a 
comparison  of  the  experimentally  v  'd  parameters  with  those  from  the  numerical  sim¬ 
ulation  at  the  final  time.  On  the  coa  ,rid,  the  total  dissipation  rate  that  was  calculated 

is  still  slightly  below  the  value  measu.i^  oy  CBC.  A  96^  grid  generates  values  of  e  consistent 
with  CBC. 


11 


CBC 

on 

9"^ 

^rmi 

5.03 

5.18 

5.19 

5.21 

|fr(Sk) 

- 

-0.42 

-0.51 

-0.52 

E 

38.6 

40.3 

40.4 

40.7 

e 

154.4 

151.3 

154.6 

156.8 

■^11 

0.34 

0.33 

0.34 

0.33 

Ai2 

- 

0.24 

0.24 

0.23 

^13 

- 

0.24 

0.24 

0.23 

Rx 

36.6 

37.8 

40.4 

37.2 

Table  2;  Final  conditions  (t=0.1145)  based  on  CBC  experiment  and 
Clark  et  al.  (1979)  calculation.  Mach  number  is  zero. 


At  t=0.1 145,  the  diagonal  components  of  Sk  are  —0.5  which  agree  well  with  the  numerical 
results  of  Kerr  (1985).  Kerr  studied  isotropic,  turbulent  flow,  but  prevented  the  decay  of 
energy  by  using  an  exterior  energy  source  at  the  large  length  scales. 

As  noted  in  the  previous  section,  we  have  cho.sen  not  to  de-alias  the  advcction  terms. 
In  reaching  this  decision  we  drew  upon  the  extensive  evidence  that  has  accumulated  on 
aliasing  effects  in  the  last  dozen  years  (Canute,  et  al  1987.,  Chapters  3,  4  and  7)  and  upon 
tests  conducted  with  the  incompressible  isotropic  turbulence  code.  In  this  code,  de-aliasing  is 
accomplished  by  applying  the  2/3  rule  (Canuto,  et  al  1987,  Chapters  3  and  7)  in  an  isotropic 
fashion;  e.g.,  the  de-aliased  results  for  a  64^  grid  are  obtained  by  running  the  incompressible 
code  on  a  96^  grid  and  applying  the  truncation  given  by  (58)  with  A7/3  in  place  of  jV/2  on 
the  right  hand  side.  The  results  are  summarized  in  Fig.  1.  Here  we  present  the  energy 
spectra  E{k)  (defined  in  the  Appendix)  for  64^,  96^,  and  128^  grids  at  t  =  0.0586  for  both 
aliased  and  de-aliased  calculations.  Some  adverse  effects  of  aliasing  are  apparent  on  the  64^ 
grid,  but  they  are  only  in  the  tail  of  the  spectra,  and  they  are  already  insignificant  on  a  96^ 
grid.  For  the  reasons  outlined  here,  a  96^  grid  was  chosen  as  the  standard  discretization  for 
the  incompressible  and  for  the  compressible  simulations. 


5.  Compressible  Turbulence  Results 

5.1.  Direct  Simulations 

Recent  work  on  the  direct  numerical  simulation  of  homogeneous  compressible  turbulence 
has  indicated  the  crucial  role  played  by  the  initial  conditions.  Passot  and  Pouquet  (1986) 
conducted  direct  simulations  of  two-dimensional,  compressible  isotropic  turbulence  and  con¬ 
cluded  that  when  the  initial  rms  density  fluctuations  are  small,  the  turbulence  statistics 
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remain  quasi-incompressible  for  turbulent  Mach  numbers  Mt  less  than  0.3. 

i)' 


=  (-) 


(68) 


They  also  demoustrated  (through  the  use  of  direct  numerical  simulations)  that  eddy  shock- 
lets  result  for  sufficiently  high  initial  rms  density  fluctuations  and/or  turbulent  Mach  num¬ 
bers.  A  more  systematic  analysis  and  categorization  of  the  effect  of  the  initial  conditions  on 
compressible  isotropic  turbulence  was  achieved  recently  by  Erlebacher  et  al.  (1990). 


They  concluded  that  for  0  <  <  0.3,  Prm»  must  initially  be  0{Mt)  for  the  resulting 

turbulence  statistics  to  become  strongly  compressible  with  an  0(1)  ratio  of  compressible 
to  incompressible  turbulent  kinetic  energy.  (For  the  range  of  Mt  considered  herein,  no 
eddy-phocklets  occur.)  On  the  other  hand,  if  initially  prmi,  Trmi  Mt,  then  the  resulting 
turbulence  statistics  remain  quasi-incompressible. 


We  first  present  the  results  of  direct  numerical  simulations  of  compressible  isotropic 
turbulence  corresponding  to  the  initial  conditions  of  the  CBC  experiment  but  with  a  variety 
of  non-zero  mean  Mach  numbers.  Since  for  these  simulations,  the  initial  conditions  are 
Prjnt  =  0,  Trm§  I)  o^ily  Weakly  comprcssible  turbulence  statistics  are  expected  according 
to  the  theoretical  results  of  Erlebacher  et  al.  (1990).  Unless  specified  otherwise,  a  subscript 
rms  for  any  variable  T  refers  to  the  quantity  I {^) ■ 

The  initial  pressure  distribution  over  the  entire  field  is  specified.  The  fluctuating  pressure, 
p/  is  determined  from  the  velocity  distribution  by  enforcing  a  zero  initial  time  derivative  for 
V  •  V.  A  Poisson  equation  for  p/  is  obtained  from  the  divergence  of  the  momentum  equation 
after  setting  the  time  variation  of  V  •  v  to  zero  (Feiereisen  et  al.  1981).  The  mean  pressure, 
Pm,  is  then  determined  so  that  a  prescribed  initial  mean  average  Mach  number,  AIq,  defined 
to  be  the  ratio  of  rms  fluid  velocity  and  rms  speed  of  sound,  is  achieved.  An  analytic 
expression  for  Pm  is  given  by 

_  Mq  Jv'^d^x  +  Jjpjp-'^d^x 
jfp-^d^x 

The  initial  average  Mach  number  is  specified  at  the  outset  of  the  direct  numerical  simulations 
(DNS)  as  an  initial  condition.  Density  is  initially  set  to  unity,  while  the  temperature,  if 
required,  is  derived  from  the  equation  of  state.  Direct  numerical  simulations  are  performed 
for  Mq  =  0.0, 0.1, 0.4  and  0.6. 

The  Mach  0.6  case  contains  localized  regions  of  supersonic  flow  as  evidenced  by  tables 
3-4  and  by  the  three-dimensional  contour  of  Mach  1  furnished  in  figure  2.  Nonetheless,  the 
statistical  properties  of  the  flow  remain  largely  unaffected  by  compressibility  effects.  This  is 
shown  in  figures  3-6  which  track  the  time  histories  of  several  statistical  variables  obtained 
from  96^  DNS.  The  time  histories  for  skewness  (fig.  3),  Aj,  (fig.  4),  and  total  kinetic  energy 
(fig.  5)  at  Mach  numbers  0.0,  0.1,  0.4  and  0.6  are  virtually  superimposed  on  each  other. 

Flatness  and  skewness  are  affected  the  most  by  compressibility  effects  in  these  simula¬ 
tions.  F’igure  3  indicates  that  the  skewness  which  corresponds  to  an  isotropic  turbulent  state 
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monotonically  increases  with  Mach  number.  It  is  -0.50  at  Mq  =  0  and  has  increased  to  -0.46 
at  Mq  =  0.6.  Before  the  flow  has  reached  a  state  of  isotropic  turbulence,  the  time  evolu¬ 
tion  of  skewness  at  all  Mach  numbers  are  indistinguishable  from  each  other.  The  physical 
system  has  equilibrated  after  approximately  one  third  the  total  computation  time.  While 
not  reaching  an  equilibrium  value,  it  is  nonetheless  worthwhile  to  point  out  that  the  flatness 
parameter  decreases  by  2%  as  the  Mach  number  is  raised  from  0.0  to  0.6  as  seen  in  figure  6. 

Figure  5  illustrates  the  decay  of  turbulent  kinetic  energy  (/  as  a  function  of 

time.  This  decay  is  a  natural  consequence  of  viscous  damping.  After  a  brief  initial  increase, 
R\  continuously  decreases  in  time,  (fig.  7),  with  no  sign  of  stabilizing.  On  the  other  hand,  \ 
which  is  representative  of  the  smaller  eddies,  increases  in  time  (fig.  4).  This  indicates  that 
energy  in  the  higher  wavenumbers  is  being  depleted  by  the  molecular  viscosity. 

Tables  3-5  summarize  the  results  of  direct  simulations  of  compressible  isotropic  turbu¬ 
lence  for  Mq  =  0.1, 0.4  and  0.6,  for  t=0.1145  on  three  different  grids^.  Incompressible  results 
are  included  for  comparison.  On  all  the  grids,  the  compressible  data  converges  to  the  in- 


Mq 

E 

e 

1^  ■ 

(itr(Sk)) 

(M) 

^max 

0.0 

40.26 

157.4 

0.00 

-0.424 

0.00 

0.00 

0.1 

40.82 

158.2 

0.17 

-0.440 

0.07 

0.21 

0.4 

41.09 

160.4 

1.50 

-0.428 

0.28 

0.84 

0.6 

41.32 

162.3 

3.30 

-0.406 

0.43 

1.26 

Table  3:  Summary  of  direct  simulations  on  a  64^  grid  with  initial  aver¬ 
age  Mach  numbers  of  0.0,  0.1,  0.4  and  0.6  at  t=0.1145. 


Mq 

E 

€ 

|V  • 

(tMSk)) 

(M) 

Mfnax 

0.0 

40.35 

154.3 

0.00 

-0.506 

0.00 

0.00 

0.1 

40.49 

155.5 

0.14 

-0.505 

0.07 

0.23 

0.4 

40.79 

157.0 

1.17 

-0.493 

0.28 

0.93 

0.6 

41.04 

158.3 

2.82 

-0.477 

0.42 

1.39 

Table  4:  Summary  of  direct  simulations  on  a  96^  grid  with  initial  aver¬ 
age  Mach  numbers  of  0.0,  0.1,  0.4  and  0.6  at  t=0.1145. 


compressible  results  as  the  Mach  number  is  driven  towards  zero.  As  expected,  the  divergence 
of  velocity  no  longer  vanishes,  and  is  now  an  increasing  function  of  Mq. 

DNS  at.  Mo  —  0.6  was  not  performed  on  the  126^  grid. 
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Table  5:  Summary  of  direct  simulations  on  a  128’  grid  with  initial 
average  Mach  numbers  of  0.0,  0.1  and  0.4  at  t=:0.1145 


While  the  dissipation  is  approximately  the  same  on  the  two  finer  grids,  the  consistently 
lower  values  on  the  coarsest  grid  confirm  the  previously  stated  conclusion  that  a  64’  grid 
cannot  resolve  all  the  length  scales.  As  a  function  of  increasing  Mach  number,  the  trace  of 
Sk  increases,  the  dissipation  decreases,  while  the  total  kinetic  energy  increases  very  slightly. 

The  results  in  tables  3-5  are  averages  over  several  DNS  runs  with  different  initial  seeds. 
A  given  seed  uniquely  determines  the  initial  velocity  distribution,  and  therefore  the  pressure 
and  temperature  fields.  Variations  of  the  seed  are  only  necessary  to  eliminate  the  statistical 
uncertainty  due  to  the  random  velocity  distribution.  The  distribution  of  velocity  on  two 
different  grid  sizes  are  different  even  when  the  initial  seed  is  the  same. 

Skev/ness  is  even  more  sensitive  to  the  grid  refinement  than  the  dissipation  as  witnessed 
by  its  decrease  from  a  value  of  -0.505  to  one  of  -0.521  on  96’  and  128’  grids  respectively.  This 
might  be  a  result  of  the  greater  sensitivity  of  the  fluctuating  velocity  field  spatial  derivatives 
to  slight  inaccuracies  in  the  flow  variables. 

5.2.  Data  Analysis 

Using  the  data  generated  from  the  previously  discussed  DNS  of  compressible  homogeneous 
turbulence  at  low  Mach  numbers,  the  proposed  subgrid-scale  (SGS)  model  is  now  validated. 
Models  relate  the  subgrid-scale  stresses  -  which  are  not  available  to  a  large-eddy  simulation 
code  -  to  the  large  scale  velocities  which  are  known.  These  velocities  are  simply  the  Favre- 
filtered  velocities  introduced  earlier.  The  Favrc-filteied  velocities  are  calculated  by  filtering 
the  resolved  DNS  velocity  field  with  a  Gaussian  spatial  filter  of  width  A  =  A/Ai/,  where 
Axf  is  the  grid  spacing  on  the  fine  grid.  For  convenience,  Ac  and  A/  refer  to  the  filter  width 
A  non-dimensionalized  with  respect  to  the  coarse  and  fine  grid  spacing  respectively. 

Perturbed  velocity  fluctuations  on  the  fine  grid  are  the  difference  between  the  fully  re¬ 
solved  velocity  and  the  filtered  ones,  given  by 

v'  =  V  -  V.  (70) 

From  v'  and  v,  subgrid-scale  stresses  based  on  DNS,  refeied  to  as  exact,  are  calculated. 
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These  include  the  Leonard,  cross  and  Reynolds  subgrid-scale  stresses  given  by  Eqs.  (20)- 
(22).  However,  these  subgrid-scale  stresses  themselves  do  not  directly  affect  the  evolution  of 
the  system.  The  momentum  equation  is  only  influenced  by  the  divergence  of  the  subgrid-scale 
stresses  (i.e.,  the  vector  level).  Similarly,  v  •  (V  •  t)  (the  scalar  level)  is  a  better  represen¬ 
tation  of  the  dissipation  terms  in  the  energy  equation  than  are  the  stresses.  Consequently, 
correlations  are  performed  on  the  tensor,  vector  and  scalar  levels.  Ideally,  high  correlations 
are  desired  on  all  levels. 

The  data  analysis  proceeds  in  multiple  stages.  First,  the  exact  stresses  calculated  from 
the  DNS  are  injected  down  to  the  coarse  grid,  along  with  the  filtered  velocities.  The  modeled 
subgrid-scale  stresses  are  then  calculated  (excluding  the  model  constants)  on  the  coarse  grid. 
Some  variables  must  be  filtered  a  second  time  (e.g.  cross  stress  terms).  Rather  than  calculate 
them  on  the  fine  grid  (which  is  not  available  to  the  large-eddy  simulation  codes),  a  Gaussian 
filter  is  applied  to  v  on  the  coarse  grid  with  a  filter  width  of  A,.  Consistency  between  the 
coarse  and  fine  filter  widths  is  achieved  by  insuring  that 

Ac  Nc 

where  Nj  and  Nc  are  respectively  the  number  of  nodes  along  one  direction  of  the  fine  and 
coarse  grids.  This  guarantees  that  the  filtering  on  the  coarse  and  fine  grids  is  performed 
over  the  same  region  in  physical  space.  Derivative  evaluations  on  both  the  coarse  and  the 
fine  grid  are  based  on  Fourier  collocation.  Calculations  by  McMillan  and  Ferziger  (1979) 
indicate  that  the  model  constants  are  sensitive  to  the  accuracy  of  the  derivative  evaluations. 
A  general  trend  that  has  been  observed  is  that  the  Smagorinsky  constant  is  lowered  when 
derivative  quantities  are  evaluated  more  accurately.  Our  constants  are  therefore  expected  to 
lie  in  the  lower  range  of  the  values  obtained  by  McMillan  (1980). 

Next,  the  model  constant,  Cr,  is  calculated.  Unfortunately,  the  constants  Ci.n  be  calcu¬ 
lated  by  a  wide  variety  of  algorithms,  each  with  its  own  merits.  Moreover,  for  each  algorithm, 
the  constants  can  be  evaluated  from  tensor,  vector  or  scalar  information.  Therefore,  criteria 
mu.st  be  established  to  identify  the  best  method.  A  key  test  is  that  L-f-C  should  be  Galilean 
invariant.  To  make  use  of  this  fact,  an  additional  constant,  Cc,  is  introduced  as  an  extra 
factor  in  the  subgrid-scale  cross  stress  model.  A  self-consistent  method  of  calculating  the 
constants  must  reproduce  Cc  =  1  to  satisfy  the  Galilean  invariance  property  stated  above 
(Speziale  1985).  Additional  tests  are  performed  on  coarse  grids  with  varying  degrees  of 
refinement  which  further  decrease  the  number  of  choices.  A  thorough  discussion  of  model 
constants  is  the  subject  of  the  next  subsection. 

Once  a  single  or  a  multiple  set  of  model  constants  have  been  determined,  the  model 
subgrid-scale  stresses  are  calculated  and  correlated  with  tl  ;  exact  subgrid-scale  stresses 
calculated  from  the  DNS  after  injection  onto  the  coarse  grid.  The  correlations  are  performed 
for  each  type  of  subgrid-scale  stress  individually,  and  for  the  total  stress  (L  -|-  C  -|-  R).  Strong 
differences  in  the  correlation  coefficients  relating  total  stresses  are  noticed  depending  on 
whether  or  not  the  Leonard  stresses  are  included.  Finally,  the  correlations  obtained  from 
the  proposed  model  are  compared  with  the  linear  combination  model,  which  has  been  shown 


16 


to  be  one  of  the  best  models  available  for  incompressible  isotropic  turbulence.  Correlation 
coefficients  are  calculated  based  on  the  two  pairs  of  constants  that  are  obtained  from  the 
above  considerations.  The  set  that  is  finally  retained  corresponds  to  the  highest  levels  of 
correlation  of  the  total  stress  on  the  vector  and  scalar  levels.  These  matters  are  treated  more 
completely  in  a  later  subsection. 

To  avoid  a  possible  confusion  of  terminology  when  referring  to  variables  being  compared 
against  each  other,  superscripts  m  and  e  are  sometimes  used.  They  respectively  refer  to 
modeled  and  exact  (based  on  DNS)  variables  at  the  tensor,  vector  and  scalar  levels. 


5.2.1.  Model  Constants 

The  proposed  model  given  by  equations  (31)  and  (42)  has  two  undetermined  coefficients. 
The  constant,  Ca,  is  associated  with  the  modeled  subgrid-scale  Reynolds  stress,  R,  w’hile 
Prx  is  associated  with  the  thermal  heat  flux. 

Although  the  cross  stress  model  has  no  constant  associated  with  it,  it  is  nonetheless 
multiplied  by  a  constant  Cc-  This  is  done  in  the  hope  of  reducing  the  number  of  schemes 
by  which  the  constants  can  be  calculated.  A  good  model  should  reproduce  a  cross  stress 
constant  of  one  to  guarantee  Galilean  invariance.  Once  the  constants  have  been  determined, 
Cc  is  set  to  one  and  forgotten.  Because  the  flow  is  isotropic,  constants  are  expected  to 
be  the  same  for  the  three  diagonal  stress  components,  the  three  off-diagonal  components 
and  the  three  vector  components.  Therefore,  the  values  presented  in  the  tables  below  are 
averaged  over  the  appropriate  components.  In  the  tables,  D  refers  to  averaged  diagonal 
components,  OD  to  averaged  off-diagonal  components,  V  to  averaged  vector  components 
and  S  to  averaged  scalar  components.  Similar  averagings  are  performed  for  the  correlation 
coefficients. 

The  two  constants  {Cn  and  Cc)  sre  calculated  using  two  techniques  -  each  applied  on 
the  tensor,  vector  and  scalar  levels.  One  method  enforces  equality  of  the  rms  levels  of  the 
exact  and  modeled  stresses.  This  is  done  for  each  individual  subgrid-scale  stress  (i.e.,  the 
subgrid-scale  Reynolds  stress  and  the  cross  stress).  Hereafter  this  approach  is  referred  to 
as  RMS.  The  second  method  utilizes  multiple  linear  least  square  regression  (LSQ)  between 
the  exact  (18)  and  the  modeled  (31)  total  subgrid-scale  stresses  to  determine  the  constants. 
Table  6  summarizes  the  results  obtained  from  incompressible  data.  The  three  cases  presented 
are  identical  except  for  the  initial  random  number  seed.  The  constants  are  independent  of 
the  detailed  velocity  statistics.  These  results  are  based  on  a  vector  level  comparison  between 
the  modeled  and  exact  stresses.  Both  RMS  and  LSQ  produce  Cc  close  to  unity  as  required. 
Unfortunately  this  prevents  a  rational  choice  from  being  made  between  the  two  approaches, 
A  more  complete  set  of  LSQ  constants  is  presented  in  table  7.  They  are  computed  at  Mach 
numbers  of  0.0,  0.1,  0.4  and  0.6  on  the  tensor,  vector  and  scalar  level.  Computations  are 
performed  on  a  coarse  grid  of  16^, 
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LSQ 

RMS 

Seed 

Cr 

Cc 

Cr 

Cc 

1 

IlMM 

■ffiMUJ 

BEa 

2 

Hg 

119 

3 

0.012 

lEi 

0.022 

1.03 

Table  6:  Model  constants  calculated  by  LSQ  and  RMS  between  exact 
and  modeled  total  stresses  (L  +  C  +  R).  Results  are  based  on  three  iden¬ 
tical  incompressible  simulations  except  for  the  random  initial  velocity 
distributions.  Calculations  are  on  the  vector  level  on  a  16^  grid. 


At  first  glance,  Cc  is  near  unity  at  both  the  scalar  and  vector  levels.  However  whereas 
on  the  vector  level,  the  constant  remains  within  4%  of  unity  for  all  Mach  numbers,  this  is 
not  the  case  on  the  scalar  level  where  Cc  is  a  decreasing  function  of  Mach  number.  This 
trend  is  present  in  the  data  generated  from  both  random  seeds.  Although  not  presented 
here,  the  rms  cross  stress  model  constant  is  also  near  unity  when  calculated  based  on  vector 
level  stresses.  Therefore,  a  preferred  method  for  the  determination  of  the  model  constants 
is  still  not  possible. 


5.2.2.  Filter  Width  and  Grid  Coarseness 

Confirmation  of  the  numerical  evidence  presented  by  McMillan  and  Ferziger  (1979)  that 
Ac  =  2  is  the  best  filter  width  is  given  in  table  8  [Mq  =  0.1).  The  criterion  used  to 
determine  the  validity  of  the  filter  width  is  that  must  remain  close  to  unity  on  the 

vector  level.  Only  when  Ac  =  2  is  Cc  near  one.  Similar  results  hold  for  LSQ  constants.  The 
constants  also  vary  with  respect  to  the  coarse  grid  on  which  the  LES  is  to  be  performed. 
Table  9  summarizes  the  model  LSQ  and  RMS  constants  evaluated  from  modeled  stresses 
calculated  on  16^  and  32^  grids  on  the  vector  level.  The  data  (Mq  =  0.1)  shows  that  Cr 
varies  by  30%  when  the  grid  size  on  which  the  modeled  stresses  are  calculated  ranges  from 
16^  to  32^.  On  the  other  hand,  large  eddy  simulations  using  finite-difference  algorithms 
might  be  performed  on  grids  as  large  as  128^.  Unless  a  subgrid-scale  model  is  found  which 
produces  constants  independent  of  the  coarse  grid  size,  LES  simulations  will  run  the  risk 
of  producing  the  wrong  results.  Perhaps  a  more  complicated  dependence  of  the  modeled 
subgrid-scale  Reynolds  stresses  on  A  is  required. 

The  best  model  constants  will  produce  the  highest  correlations  between  the  modeled  and 
exact  subgrid-scale  stresses  on  all  levels  (tensor,  vector  and  scalar).  Based  on  the  previous 
discussion,  only  constants  calculated  on  a  vector  level  are  adequate  because  they  produce 
a  Cc  of  unity.  Unfortunately,  a  clear  cut  choice  between  RMS  and  LSQ  constants  cannot 
be  made  because  Cc  (on  the  vector  level)  is  nearly  one  in  both  cases.  Rather  than  make 
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an  arbitrary  choice,  both  sets  of  constants  are  considered  when  correlating  exact  against 
modeled  stresses.  For  reference,  the  constants  used  henceforth  are 


=  0.012  (72) 

=  0.023.  (73) 


Seed  1 

Seed  2 

Mo 

0.0 

0.1 

0.4 

0.6 

0.0 

0.6 

D 

1.32 

1.32 

1.32 

1.31 

1.32 

1.31 

Cc 

V 

1.04 

1.04 

1.02 

1.00 

1.04 

1.00 

S 

1.00 

1.00 

0.95 

0.873 

0.971 

0.934 

D 

0.018 

0.018 

0.018 

0.018 

0.016 

0.015 

Ch 

V 

0.012 

0.012 

0.012 

0.012 

0.012 

0.012 

S 

0.015 

0.015 

0.015 

0.014 

0.015 

0.015 

Table  7;  LSQ  model  constants.  Filter  widths  are  Ay  —  12  and  Ac  =  2. 


LSQ 

RMS 

A/ 

6 

12 

24 

6 

12 

24 

Ac 

1 

2 

4 

1 

2 

4 

Ch 

0.007 

0.012 

0.020 

0.019 

0.023 

0.034 

Cc 

0.31 

1.03 

1.33 

0.82 

1.03 

1.13 

Table  8:  LSQ  and  RMS  model  coefficients  between  exact  and  modeled 
stresses  on  the  vector  level  at  Mo  —  0.1.  Results  are  obtained  with  fine 
filter  widths  ot  u,  12  and  24  while  maintaining  the  proper  ratio  of  6 
between  fine  and  coarse  widths.  The  coarse  grid  is  16^. 


5.2.3.  Correlations 

Correlation  coefficients  have  long  been  a  preferred  diagnostic  tool  for  estimating  the  reliability 
of  the  modeled  stresses.  However,  the  subgrid-scale  stresses  have  been  separated  into  various 
components  (L,C,R),  each  modeled  separately  and  although  the  correlation  of  any  of  these 
components  against  their  models  might  be  excellent,  it  is  still  possible  for  the  correlation 
of  the  exact  total  stiess  against  the  modeled  total  stress  to  be  less  impressive.  Such  is  the 
case,  for  instance,  when  two  stress  components  have  opposite  signs  and  partially  cancel  each 


19 


LSQ 

RMS 

Grid 

Cr 

Cc 

Cr 

Cc 

16" 

1.03 

32" 

1.03 

in 

Table  9:  LSQ  and  RMS  model  constants  based  on  subgrid-scale  stresses 
evaluated  on  16^  and  32^  coarse  grids. 


other  out.  As  a  final  note,  before  the  specific  correlation  coefficients  are  presented,  one  must 
always  be  attentive  to  the  actual  relationship  between  the  exact  and  the  modeled  variable, 
even  when  the  correlation  coefficient  is  relatively  high  (say,  above  70%).  A  correlation 
coefficient  as  high  as  70%  may  not  be  as  good  as  it  seems.  For  example,  the  correlation 
between  the  functions  y  =  x  and  y  =  exp(— i)  in  the  interval  [0, 1]  is  approximately  70% 
although  they  are  qualitatively  different  functions!  As  a  consequence,  correlations  are  only 
deemed  good  if  the  correlation  coefficient  is  above,  say,  the  90%  level. 

For  convenience,  the  compressible  subgrid-scale  model  is  restated  here; 


oR.,  =  2CflpA2  CSkiSkif"' 

(74) 

c. 

II 

o 

(76) 

Cij  -  -p{vtVj  -  v.Vj) 

(76) 

The  correlation  coefficients  presented  in  table  10  between  exact  and  modeled  R  and  C, 
are  independent  of  the  model  constants.  The  correlation  coefficients  are  insensitive  to  the 
average  Mach  number  variation.  The  neglect  of  appears  to  be  a  good  approximation;  the 
direct  simulations  .show  it  to  be  several  orders  of  magnitude  smaller  than  the  thermodynamic 
pressure.  For  example,  for  all  of  the  compressible  isotropic  turbulence  simulations  conducted 
in  this  study, 

^  <  3  X  10-" 

rm# 

and,  hence,  the  effect  on  the  isotropic  part  of  the  Reynolds  stress  tensor  is  dominated  by  the 
thermodynamic  pressure. 

In  most  of  the  literature  on  subgrid-scale  models,  the  Leonard  stress  has  been  omitted 
from  the  total  stress  correlations  on  the  grounds  that  it  is  calculated  exactly  (Bardina, 
Ferziger,  and  Reynolds  1983,  McMillan  1980,  McMillan  and  Ferziger  1979).  However,  it  has 
recently  been  shown  by  Speziale  (1985)  that  the  combination  L  -f  C  is  Galilean  invariant, 
while  the  Leonard  and  cross  stresses,  individually,  are  not.  Therefore  we  feel  that  correlations 
of  the  total  stress  should  include  the  Leonard  stress.  Table  11  summarizes  the  correlation 
coefficients  between  the  total  stress  including  and  excluding  the  Leonard  stress.  Results 
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Mo 

0.1 

0.4 

0.6 

El 

31 

31 

31 

R 

lail 

26 

26 

25 

11 

22 

22 

22 

n 

45 

45 

45 

D 

89 

89 

89 

C 

OD 

91 

91 

91 

V 

80 

80 

80 

s 

75 

74 

72 

Table  10:  Correlations  between  exact  and  modeled  stresses,  R  and  C, 
at  Mach  0.1,  0.4  and  0.6. 


are  presented  at  Mach  0  and  Mach  0.6.  When  the  Leonard  stresses  are  left  out,  correlation 
coefficients  similar  to  those  of  Bardina  are  obtained  on  all  levels.  However,  the  inclusion  of  L 
decreases  the  correlations  at  the  vector  level  by  approximately  30%.  Table  11  substantiates 
that  the  correlation  coefficients  (with  and  without  L)  are  nearly  independent  of  the  initial 
average  Mach  number. 

Correlation  coefficients  between  (C  +  R)*  and  various  combinations  of  the  modeled 
stresses  using  constants  calculated  by  LSQ  and  RMS  are  summarized  in  table  12.  The 
second  column  indicates  the  model  against  which  the  total  stress  is  being  compared.  Al¬ 
though  at  first  glance  RMS  based  constants  perform  better  at  the  tensor  level  when  the 
total  modeled  stress  is  considered,  at  the  vector  and  scalar  levels,  the  trend  is  reversed. 
Correlations  at  the  vector  and  scalar  level  are  higher  by  4%  using  the  LSQ  constants. 

When  the  constants  are  selected  based  on  LSQ,  table  12  indicates  that  the  correlations 
on  all  levels  are  highest  when  all  modeled  components  are  included.  However,  the  dynamic 
evolution  of  the  large  scale  velocities  only  brings  into  play  the  stresses  on  the  vector  and 
scalar  level.  Therefore  the  coefficients  which  produce  the  maximum  correlations  of  total 
stress  on  these  two  levels  should  be  chosen.  This  leads  to  an  optimum  choice  of 

Ch  =  0.012  (78) 

calculated  by  least  squares  fit  of  the  total  stress  on  the  vector  level. 

Conversion  of  Cr  to  the  standard  form  currently  used  in  incompressible  LES  ^  produces 
a  Smagorinsky  constant  of  0.092.  McMillan  (1979)  obtained  a  value  of  Cs  =  0.10  when 
spectral  collocation  derivative  computations  were  employed.  This  constant  corresponds  to 
a  scalar  level  evaluation  based  on  RMS.  On  the  vector  level,  McMillan  calculated  a  higher 

^In  a  number  of  reports  (McMillan  and  Ferzigcr  1979,  Bardina  et  al.  1983),  the  subgrid-scale  Reynolds 
stress  model  is  proportional  to  C|.  In  these  cases,  the  relationship  between  Cs  and  Cjf  is  C/j  =  v/2Cj. 
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Mo  = 

0.0 

Mo  — 

0.6 

L+C+R 

C-fR 

L  4-  C  "h  R 

C+R 

D 

93 

82 

93 

81 

OD 

80 

85 

79 

84 

46 

72 

46 

71 

5 

56 

73 

56 

74 

Table  11:  Comparison  of  correlation  coefficients  of  the  exact  total  stress 
versus  its  model.  The  modeled  terms  are  computed  on  a  16^  coarse  grid. 
Each  case  is  presented  with  and  without  the  inclusion  of  the  Leonard 
subgrid-scale  stress  terms  (calculated  exactly).  Both  the  incompressible 
and  the  M=0.6  case  are  shown  to  illustrate  the  weak  influence  of  Mach 
number  on  the  correlation  coefficients  with  Cr  —  0.012. 


Least  squares 

RMS 

Cr  -  0.023 

Exact 

Model 

D  OD  V  S 

D  OD  V  S 

C-fR 

C 

C  +  R 

78  82  68  62 

81  84  71  70 

78  82  68  62 
88  83  67  67 

Table  12:  Correlation  of  the  exact  total  stress  (C  -f  R)*  with  various 
models.  The  modeled  stresses  are  defined  in  equations  (74)-(76). 


Smagorinsky  constant  of  0.13.  This  value  can  be  obtained  from  the  present  data  by  using 
instead  of  oR^^^-  However  as  shown  above,  the  correlations  of  the  total  subgrid- 
scale  stress  would  be  lower. 

Initial  tests  of  the  subgrid-scale  heat  flux  model  produced  a  turbulent  Prandtl  number 
in  the  range  of  0.4  to  0.5.  For  this  initial  study,  we  take 

Ptt  =  0.5  (79) 

which  is  a  value  that  has  been  used  in  previously  published  large-eddy  simulations  of  turbu¬ 
lent  flows  with  thermal  convection  (cf.  Eidsor  1985),  A  more  accurate  calculation  of  Prj 
could  be  accomplished  in  a  compressible  flow  with  significant  mean  temperature  gradients; 
this  is  beyo  id  the  scope  of  the  present  study  and  must  await  future  research. 
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6.  Large-Eddy  Simulation  of  C':.*ipressible  Isotropic  Turbulence 


Now,  in  order  to  demonstrate  the  efficacy  of  the  subgrid-scale  models  derived  in  this  paper, 
a  large-eddy  simulation  of  compressible  isotropic  turbulence  is  conducted.  Since  most  high¬ 
speed  compressible  turbulent  flows  have  significant  regions  where  the  turbulence  statistics  are 
quasi-incompressible,  it  is  important  that  the  models  perform  well  for  weakly  compressible 
turbulence  -  the  type  of  flow  considered  in  the  last  section  for  the  a  priori  tests.  However, 
it  is  well  known  that  a  priori  tests  only  provide  a  relatively  weak  gauge  for  the  performance 
of  subgrid-scale  models  in  an  actual  large-eddy  simulation  (see  Hussaini,  Speziale  and  Zang 
1990).  It  is  therefore  important  to  examine  their  performance  in  an  actual  large-eddy  sim¬ 
ulation,  particularly  for  a  case  where  significant  compressibility  eflfects  are  exhibited  in  the 
turbulence  statistics. 

Direct  simulations  of  compressible  isotropic  turbulence  were  conducted  corresponding  to 
the  initial  conditions  Re  =  250,  =  0.1,  (prm.)o  =  0,  {Trn$)o  =  0.0626  and  two  values 

of  Xo^  0  and  0.2  which,  respectively,  correspond  to  initial  values  for  {Rx)o  of  30.0  and  26.3. 
Here,  x  =  /{E^  -I-  E^)  where  E'  and  E’^  are  the  incompressible  and  compressible  parts  of 

the  turbulent  kinetic  energy,  respectively.  The  direct  numerical  simulations  were  performed 
on  a  96^  mesh.  Characteristic  energy  and  dissipation  spectra  associated  with  the  DNS  are 
shown  in  figures  8  (a)-(c)  and  9  (a)-(c)  for  Xo  —  0  and  xo  ~  0-2,  respectively.  Figure  8  (b) 
clearly  shows  that  x  remains  very  small  for  all  time,  with  only  very  slight  modification  of  the 
energy  spectrum  in  time.  In  contrast,  figure  9  (b)  shows  an  initial  cascade  of  the  spectrum 
towards  the  smaller  scales,  followed  by  a  strong  energy  dissipation  at  the  smaller  scales  of  up. 
The  dissipation  spectra  {k'^E{k))  shown  in  figures  8  (a)  and  9  (a)  demonstrate  that  both  the 
small  and  large  scales  are  well  resolved.  Comparing  figures  8  (c)  and  9(c),  the  incompressible 
energy  E^  is  the  same  at  t  =  0  and  t  =  1.6,  but  is  influenced  by  compressibility  effects  at 
the  intermediate  times.  This  influence  is  characterized  by  a  slight  decrease  in  E^  for  Xo  >  0. 

Integral  properties  of  these  two  cases  are  plotted  in  figures  10  (a)-(e).  These  figures 
contrast  the  two  runs  corresponding  to  Xo  =  0  and  Xo  =  0.2.  The  higher  compressibility 
has  a  variety  of  effects  on  the  flow.  The  total  kinetic  energy  decays  at  a  slightly  slower  rate, 
while  both  the  skewness  and  the  flatness  are  decreased.  In  other  words,  finite  compressibility 
drives  the  flow  more  towards  a  Gaussian  state.  Both  the  integral  scale  Lj  and  the  Taylor 
microscale  An  become  smaller  as  xo  increases.  Lastly,  the  decay  rate  of  the  microscale 
Reynolds  number  is  slower  for  finite  Xo-  A  more  detailed  study  of  these  effects  awaits  future 
research.  It  would  be  particularly  useful  to  compute  these  statistical  quantities  based  on  the 
solenoidal  and  irrotational  components  of  velocity. 

As  noted  by  Frisch  and  Orszag  (1990),  the  three  dimensional  vorticity  is  formed  of 
tubular-like  structures.  This  is  graphically  represented  in  figure  11  which  shows  a  volume 
rendering  of  Volumetric  rendering  is  a  visualization  technique  whereby  rays  are  projected 
from  the  eye  through  the  flow  (which  emits  and  absorbs  light).  The  tubular  structures  of 
vorticity  are  contrasted  wfith  the  spherical  structure  of  (V  ■  v)^  shown  in  figure  12.  This 
difference  is  to  be  expected  since  the  dilatation  satisfies  an  isotropic  wave  equation  to  lead- 
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ing  order  which  shows  no  preferential  direction.  On  the  other  hand,  the  vorticity  stretching 
occurs  along  the  axis  of  the  vorticity  vector,  which  then  becomes  a  preferred  direction. 


The  results  for  the  direct  simulation  with  the  strongest  compressibility  effects  (xo  =  0.2) 
were  filtered  and  injected  onto  a  coarse  32^  grid  for  comparison  with  an  LES  which  was  also 
performed  on  a  32^  grid  using  the  subgrld-scale  models  derived  herein.  In  this  manner,  a 
direct  comparison  can  be  made  between  the  results  of  the  LES  and  the  direct  simulation 
for  a  flow  where  the  turbulence  statistics  exhibit  significant  compressibility.  The  following 
turbulence  statistics  were  compared: 

(dv  \  ^ 

j 


(2)  the  integrated  average  total  and  isotropic  vortex  stretching  (denoted  by  u),S,j(jJj  and 

•  v,  respectively), 

(3)  the  mean  turbulence  Mach  number  defined  by 

(4)  the  mean  compressible,  incompressible  and  total  turbulent  kinetic  energy  denoted  by 

and  Etot,  respectively, 

(5)  the  level  of  compressibility  x  defined  by  x  =  I  Etot,  ‘‘■nd 

(6)  the  rms  of  the  thermodynamic  pressure,  density,  and  temperature  denoted  by  Prm», 
Prmt,  and  Trm„  respectively. 

These  quantities  represent  a  good  choice  of  turbulence  statistics  to  monitor  the  effects  of 
compressibility. 

In  figure  13  (a)-(f),  a  direct  comparison  of  these  statistics  for  the  DNS  and  LES  is  made 
for  a  filter  width  A/  =  2.  It  is  clear  that  the  LES  does  an  excellent  job  in  reproducing 
the  results  of  the  DNS  with  the  possible  exception  of  x-  If  should  be  noted  that  x  exhibits 
acoustic  oscillations  and  hence  is  a  difficult  quantity  to  predict  accurately;  nonetheless,  the 
LES  yields  results  that  are  in  good  qualitative  agreement  with  the  direct  simulation.  The 
most  striking  result  is  how  well  the  compressible  turbulent  kinetic  energy  and  dilatational 
terms  are  captured. 

It  was  found  that  a  change  in  the  filter  width  A/  to  1  or  3  -  and  an  adjustment  of  the 
SGS  model  constants  Cr  and  Prr  of  up  to  25%  -  only  led  to  a  mild  degradation  c  l  the 
results.  However,  the  subgrid-scale  models  did  play  a  crucial  role  in  obtaining  the  accurate 
results  shown  in  figure  13  (i.e.,  a  32^  direct  simulation  is  substantially  under-resclved).  To 
illustrate  this  point,  the  results  of  a  direct  simulation  on  the  coarse  32^  grid  (i.e.,  an  LES 
for  A/  =  0)  is  shown  in  figure  14  for  the  same  test  case.  It  is  clear  from  this  figure  that  the 
coarse  grid  DNS  does  a  poor  job  in  capturing  the  incompressible  as  well  as  the  compressible 
turbulence  statistics.  Consequently,  it  is  the  adequate  performance  of  the  subgrid-scale 
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models  in  draining  the  proper  amount  of  energy  from  the  filtered  fields  (to  account  for  the 
presence  of  unresolved  scales)  which  leads  to  the  excellent  results  obtained  in  figure  13. 
Considering  the  degree  of  compressibility  exhibited  by  the  turbulence  statistics  for  the  test 
case  considered  in  this  section,  it  would  appear  that  the  feasibility  of  the  proposed  subgrid- 
scale  models  has  been  established.  A  more  extensive  parametric  study  of  this  subgrid-scale 
model  is  under  investigation  by  Dahlburg  et  al.  (1990). 


7.  Conclusion 


New  subgrid-scale  models  for  compressible  turbulence  have  been  developed  and  tested  against 
the  results  of  direct  numerical  simulations  of  compressible  isotropic  turbulence.  These  com¬ 
pressible  subgrid-scale  models,  which  were  based  on  the  Favre-filtered  equations  of  motion 
for  an  ideal  gas,  contain  two  dimensionless  constants  and  reduce  to  the  linear  combination 
model  of  Bardina,  Ferziger  and  Reynolds  (1983)  in  the  incompressible,  isothermal  limit.  The 
subgrid-scale  stress  model  constant  was  found  to  assume  the  value  of  Cr  —  0.012  which  gives 
rise  to  correlations  between  the  exact  and  modeled  stresses  that  were  above  70%  on  the  ten¬ 
sor,  vector  and  scalar  levels  -  a  correlation  which  compares  favorably  with  those  obtained 
in  earlier  work  on  the  subgrid-scale  modeling  of  incompressible  turbulent  flows.  Another  en¬ 
couraging  feature  lies  in  the  fact  that  these  constants  and  their  associated  correlations  were 
found  to  be  comparatively  insensitive  to  a  mean  Mach  number  in  the  range  0  <  Me  <  0.6. 

The  results  of  a  coarse  grid  32^  LES  of  compressible  isotropic  turbulence  conducted  with 
the  subgrid-scale  models  derived  in  this  paper  were  shown  to  be  in  excellent  agreement  with 
those  obtained  from  a  96'^  direct  simulation.  These  results  are  extremely  encouraging  since, 
for  the  case  considered,  on  average  25%  of  the  turbulent  kinetic  energy  was  cc  .npressible. 
Furthermore,  the  ability  for  the  LES  to  accurately  capture  the  dilatational  statistics  of  the 
flow  was  quite  surprising. 

Future  research  will  be  directed  on  several  fronts.  The  large-eddy  simulation  of  a  com¬ 
pressible  turbulent  flow  with  mean  temperature  gradients  could  lead  to  refinements  in  the 
subgrid-scale  heat  flux  model.  Furthermore,  the  large-eddy  simulation  of  compressible,  ho¬ 
mogeneous  shear  flow  could  yield  new  insights  into  the  performance  of  these  subgrid-scale 
models.  Near-wall  modifications  will  also  be  implemented  that  allow  for  the  large-eddy  sim¬ 
ulation  of  compressible,  wall-bounded  turbulent  flows.  While  further  improvements  are  still 
possible,  we  believe  that  the  essential  foundation  for  the  large-eddy  '•imulation  of  compress¬ 
ible  turbulent  flows  has  been  established  in  this  study.  With  future  research,  compressible 
LES  could  have  a  profound  impact  on  the  analysis  of  supersonic  and  hypersonic  flows  of 
aerodynamic  importance. 
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Appendix 


Initial  energy  spectrum 


Both  the  incompressible  and  compressible  direct  simulations  impose  a  specified  energy 
spectrum  on  the  initial  random  velocity  distribution.  For  comparative  purposes,  the  data 
was  obtained  from  tabular  data  found  in  Comte-Bellot  and  Corrsin  (1971).  They  tabulate 
the  function  Eii{k)  which  is  related  to  the  energy  spectrum  E{k)  by 


E{k)  = 


-k^— 
2  dk 


1  dEn 
k  dk 


(Al) 


Unfortunately  the  data  is  noisy,  so  a  least  squares  fit  is  performed  on  log£?ii,  expressed  as 
a  fourth  order  polynomial  in  log  A:.  The  final  form  obtained  for  is 


log(JE;ii )  =  2.64359  -  0.72602(log  k)  -  0.32585(log  k)^ 

+0.03525(log  ky  -  0.02344(log  k)\  ( A2) 


Calculation  of  model  constants 


Before  correlating  the  total  exact  subgrid-scale  stress  against  its  model,  the  model  con¬ 
stants  must  be  determined.  There  are  many  ways  of  accomplishing  this  among  which  two 
are  retained.  The  total  modeled  stress  is  written  as  a  linear  combination  of  modeled  terms 


T”’  =  f;c.Tr 

1=1 

while  the  exact  total  stress  is  simply 

(A3) 

(A4) 

The  unknown  constants  to  be  determined  are  the  C,.  The  first  approach  adopted  is  to 
calculate  the  root  mean  square  of  the  pairs  C,t/"  and  t*  and  equate  them  for  each  value  of 
i.  The  constants  thus  take  the  values 


(A5) 


This  method  is  referred  to  by  RMS  in  the  text. 


A  least  square  method  applied  to  the  total  stress  as  a  whole  is  an  alternative  approach. 
In  this  case,  the  norm 

I|T”-Tl|’  =  |ix:(r,'-C,T')|!’  (A6) 

1=1 

is  minimized  with  respect  to  the  coefficients  C,.  This  gives  rise  to  a  linear  system  in  the 
coefficients  which  can  be  solved  by  direct  methods  if  the  number  of  constants  is  not  too 
large.  For  the  subgrid-scale  model  considered  in  the  text,  n  =  3. 
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Definitions 

For  reference  purposes,  several  statistical  definitions  are  provided  here.  All  variables  are 
defined  on  a  three-dimensional  grid  and  are  subscripted  by  a  single  index  i  for  convenience. 
The  average  of  a  function  Ti  is 

W  =  (A7) 

As  a  function  of  the  average,  the  rms  of  !F  is 

(A8) 

Correlation  coefficients  are  fundamental  in  evaluating  subgrid-scale  models.  The  correlation 
coefficient  between  two  functions  .^and  Q  is 

C(T,Q)  =  (A9) 
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on  a  96^  grid.  CBC  initial  conditions. 

Fig.  8  Energy  and  dissipation  spectra  of  96^  DNS  at  t  =  0,0.2, 0.4, 0.8, 1.6.  Aiq  =  10,  Xo  “  0.0, 
(Mt)o  =  0.1,  {Rx)o  =  30.0,  Re  =  250. 

(a)  Dissipation,  (b)  Compressible  energy,  (c)  Solenoidal  energy. 

Fig.  9  Energy  and  dissipation  spectra  of  96^  DNS  at  f  =  0, 0.2, 0.4, 0.8, 1.6.  A:o  =  10,  xo  =  0.2, 
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(a)  Dissipation,  (b)  Compressible  energy,  (c)  Solenoidal  energy. 

Fig.  10  Turbulence  statistics  for  96^  DNS  at  <  =  0.8.  Parameters  are  Aiq  =  10,  Xo  =  0.0  and 
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I'^ig.  11  Volumetric  representation  of  Parameters  are  identical  to  those  of  Fig.  10. 

Fig.  12  Volumetric  representation  of  (V  ■  v)*.  Parameters  are  identical  to  those  of  Fig.  10. 

Fig.  13  Comparison  between  LES  results  ( - , - , - )  and  results  from  a  96^  DNS 

(A,  0,0)  injected  onto  the  32^  LES  grid.  Filter  width  Ac  =  A/  -  2. 

Fig. 14  Comparison  between  LES  results  ( - , - , - )  and  results  from  a  96^  DNS 

(A,  □,  O)  injected  onto  the  32^  LES  grid.  Filter  width  Ac  =  A/  =  0. 
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Fig.  3  Time  history  of  one  third  the  trace  of  Sk  for  Mq  of  0.0,  0.1,  0.4  and  0.6.  Direct 
simulation  was  performed  on  a  96^  grid.  CBC  initial  conditions. 
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Fig.  6  Time  history  of  one  third  the  trace  of  FI  for  Mq  of  0.0,  0.1,  0.4  and  0.6. 
simulation  was  performed  on  a  96^  grid.  CBC  initial  conditions. 
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Fig-  7  Time  history  of  Rx  for  Mq  of  0.0,  0.1,  0.4  and  0.6.  Direct  simulation  was  performed 
on  a  96^  grid.  CBC  initial  conditions. 
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Fig.  9  Energy  and  dissipation  spectra  of  96^  DNS  at  t  -  0,0.2, 0.4, 0,8, 1.6.  /cq 

(M<)o  -  0.1,  («a)o  =  26.4,  250. 

(a)  Dissipation,  (b)  Compressible  energy,  (c)  Solenoidal  energy. 
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Fig.  10  Turbulence  statistics  for  96^  DNS  at  t  =  C.8.  Parameters  are  «o  =  10,  xo  =  O.C  and 
Xo  =  0.2,  (M,)o  =  0.1,  {Prm.)o  =  0,  (T,„,)o  =  0.0526,  Re  =  250. 

(a)  Etot  =  E'^  +  FJ,  (b)  trace  of  skewness  tensor,  (c)  trace  of  flatness  tensor,  (d)  trace 
of  Taylor  length  tensor,  (e)  integral  length  scale,  (f)  microscale  Reynolds  number. 


■  ig.  11  Volumetric  representation  ofcj^.  Parameters  are  identical  to  those  of  Fig.  10. 


ig.  12  Volumetric  representation  of  (V  •  v)'^-  Parameters  are  identical  to  tho.'-.c  of  Fig.  10. 
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